

rm(list=ls(all=TRUE))

library(readstata13)
library(MASS)
library(MCMCpack)
library(xtable)



### FIGURE 1 

set.seed(123)

country <- "Germany"

load("~/Desktop/replication/Germany.rda")




quartz(type="pdf", file="~/Desktop/replication/Germany.pdf", width=6*1.62, height=6)
par(mar = c(3,3,1,1), mgp=c(2,0.5,0), family="CMU Serif")

plot(jitter(data$year), data$coopscore.postmean, type="n", cex=0.8, col="gray", xlab="Year", ylab="Relationship Score", xaxt="n")
axis(1, at=2001:2009, labels=2001:2009)
points(jitter(data$year[(data$type.1==1 & data$type.2!=1) | (data$type.1!=1 & data$type.2==1)]-0.1, factor=0.3), data$coopscore.postmean[(data$type.1==1 & data$type.2!=1) | (data$type.1!=1 & data$type.2==1)], pch=16, cex=1, col=adjustcolor("gray40", alpha.f=0.5))
points(jitter(data$year[(data$type.1!=1 & data$type.2!=1)] + 0.1, factor=0.3), data$coopscore.postmean[(data$type.1!=1 & data$type.2!=1)], pch=16, cex=1, col=adjustcolor("gray70", alpha.f=0.5))
points(data$year[data$type.1==1 & data$type.2==1], data$coopscore.postmean[data$type.1==1 & data$type.2==1], pch=16, cex=1,)

# CDU/SPD
points(data$year[(data$name.1=="CDU/CSU" & data$name.2=="SPD") | (data$name.2=="CDU/CSU" & data$name.1=="SPD")], data$coopscore.postmean[(data$name.1=="CDU/CSU" & data$name.2=="SPD") | (data$name.2=="CDU/CSU" & data$name.1=="SPD")], type="l", lwd=3, lty=1)
# Greens/SPD
points(data$year[(data$name.1=="Grune" & data$name.2=="SPD") | (data$name.2=="Grune" & data$name.1=="SPD")], data$coopscore.postmean[(data$name.1=="Grune" & data$name.2=="SPD") | (data$name.2=="Grune" & data$name.1=="SPD")], type="l", lwd=3, lty=2)
# CDU/FDP
points(data$year[(data$name.1=="CDU/CSU" & data$name.2=="FDP") | (data$name.2=="CDU/CSU" & data$name.1=="FDP")], data$coopscore.postmean[(data$name.1=="CDU/CSU" & data$name.2=="FDP") | (data$name.2=="CDU/CSU" & data$name.1=="FDP")], type="l", lwd=3, lty=3)

legend("bottomleft", c("CDU/CSU - SPD", "SPD - Grüne", "CDU/CSU - FDP"), lty=c(1, 2, 3), lwd=3)

dev.off()




## TABLE 1
# Get coefficients from stata estimation, simulate 500 draws from each, combine

rm(list=ls(all=TRUE))

data <- read.dta13("~/Desktop/replication/m1.dta")

coef <- rvec <-  NULL
for(i in 1:1000){
	usecoef <- as.matrix(data[1,(((i-1)*23)+1):(((i-1)*23)+11)])
	usevcov <- as.matrix(data[,(((i-1)*23)+12):(((i-1)*23)+22)])
	user2 <- as.matrix(data[1,(((i-1)*23)+23):(((i-1)*23)+23)])

	draw <- mvrnorm(500, usecoef, usevcov)

	coef <- rbind(coef, draw)
	rvec <- rbind(rvec, user2)
}

m1 <- as.mcmc(coef)
msum <- summary(m1)

tab <- NULL
for(i in 1:11){
	add <- matrix(c(rownames(msum$statistics)[i], format(round(msum$statistics[i,1], 3), nsmall=3, scientific=10), "", paste0("(", format(round(msum$quantiles[i,1], 3), nsmall=3, scientific=10), ", ", format(round(msum$quantiles[i,5], 3), nsmall=3, scientific=10), ")")), nrow=2, byrow=T)
	tab <- rbind(tab, add)
}

m1tab <- tab
m1tab <- rbind(m1tab, c("", round(mean(rvec), 3)))




data1 <- read.dta13("~/Desktop/replication/m2_1.dta")
data2 <- read.dta13("~/Desktop/replication/m2_2.dta")
data <- data.frame(data1, data2)

coef <- rvec <- NULL
for(i in 1:1000){
	usecoef <- as.matrix(data[1,(((i-1)*49)+1):(((i-1)*49)+24)])
	usevcov <- as.matrix(data[,(((i-1)*49)+25):(((i-1)*49)+48)])
	user2 <- as.matrix(data[1,(((i-1)*49)+49):(((i-1)*49)+49)])

	draw <- mvrnorm(500, usecoef, usevcov)

	coef <- rbind(coef, draw)
	rvec <- rbind(rvec, user2)
}

m2 <- as.mcmc(coef)
msum <- summary(m2)

tab <- NULL
for(i in 1:24){
	add <- matrix(c(rownames(msum$statistics)[i], format(round(msum$statistics[i,1], 3), nsmall=3, scientific=10), "", paste0("(", format(round(msum$quantiles[i,1], 3), nsmall=3, scientific=10), ", ", format(round(msum$quantiles[i,5], 3), nsmall=3, scientific=10), ")")), nrow=2, byrow=T)
	tab <- rbind(tab, add)
}

m2tab <- tab
m2tab <- m2tab[c(1:20, 47, 48),]
m2tab <- rbind(m2tab, c("", round(mean(rvec), 3)))




data <- read.dta13("~/Desktop/replication/m3.dta")

coef <- rvec <- NULL
for(i in 1:1000){
	usecoef <- as.matrix(data[1,(((i-1)*13)+1):(((i-1)*13)+6)])
	usevcov <- as.matrix(data[,(((i-1)*13)+7):(((i-1)*13)+12)])
	user2 <- as.matrix(data[1,(((i-1)*13)+13):(((i-1)*13)+13)])

	draw <- mvrnorm(500, usecoef, usevcov)

	coef <- rbind(coef, draw)
	rvec <- rbind(rvec, user2)	
}

m3 <- as.mcmc(coef)
msum <- summary(m3)

tab <- NULL
for(i in 1:6){
	add <- matrix(c(rownames(msum$statistics)[i], format(round(msum$statistics[i,1], 3), nsmall=3, scientific=10), "", paste0("(", format(round(msum$quantiles[i,1], 3), nsmall=3, scientific=10), ", ", format(round(msum$quantiles[i,5], 3), nsmall=3, scientific=10), ")")), nrow=2, byrow=T)
	tab <- rbind(tab, add)
}

m3tab <- tab
m3tab <- rbind(m3tab[1:10,], c("", ""), c("", ""), c("", ""), c("", ""), c("", ""), c("", ""), c("", ""), c("", ""), c("", ""), c("", ""), m3tab[11:12,])
m3tab <- rbind(m3tab, c("", round(mean(rvec), 3)))


vars <- c("Coalition", "", "Opposition", "", "Delta CMP Left-Right", "", "Delta CMP Nationalism/Multiculturalism", "", "Delta Vote Share", "",  "Election Year", "", "Number of Parties", "", "Population (log)", "", "Number of Events", "", "Mean Cooperation Score", "", "Constant", "", "R2")
outtable <- xtable(cbind(vars, m1tab[,2], m2tab[,2], m3tab[,2]))
outtable



